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Abstract. We suggest a method to experimentally obtain two-dimensional matter- 
wave discrete solitons with a self- repulsive BEC in optical lattices. At the edge of 
the Brillouin zone, a wave packet effective mass is negative which could be treated 
as inversion of the nonlincarity sign. Above critical nonlinearity this makes the wave 
packets collapse partially into localized modes with a chemical potential located in 
the gap between the first and the second bands. This critical nonlinearity is also 
associated with the smallest nonlinearity for which the discrete solitons are possible 
in the gap. Extensive numerical simulations for square and asymmetric honeycomb 
lattices in continuous model illustrate every stage of the process. 
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1. Introduction 

Periodic lattices with substantial nonlinearities appear in various systems such as 
biological molecules P, nonlinear optical wave guides [2], solid-state materials 0111, 
and Bose- Einstein condensates (BEC) jH]- In these systems, interplay between linear 
coupling effects among adjacent sites and nonlinearity can result in a self-localized state 
- lattice or 'discrete' soliton (DS) [H 13 El IH E] • Until recently direct observation of 
DS has been done only in one-dimensional (ID) optical wave guides |21 El HI- Yet 
in systems with dimensionality more than one a number of fundamental phenomena, 
such as vortex lattice solitons, bright lattice solitons that carry angular momentum, are 
expected jH] . Recently a novel experimental technique to produce photonic crystal with 
optical induction allowed the authors of Ref. P to directly observe two-dimensional 
(2D) DS. 

Dynamics of the optical pulses in nonlinear photonic crystals and BEC in optical 
lattices are governed by a nonlinear Shrodinger equation (NLSE) with periodic potential, 
hence many predictions done with respect to photonics are expected with BEC. In 
the case of BEC, the nonlinear coefficient can be either positive or negative for 
repulsive or attractive atomic interaction respectively, with most of the experiment 
being done with repulsive atoms. One-dimensional solitons in the absence of periodic 
potential were observed in attractive BECs |TIHlllj. One-dimensional matter- wave DS in 
optical lattices were studied extensively theoretically both for attractive and repulsive 
interactions ^21 ^1 Decoherence of the repulsive BEC during Bloch oscillation 
in ID optical lattices observed in experiments [T^ IT^ [T7j was related theoretically to 
generation of DS ^] (similar decoherence phenomena in two- and three-dimensional 
(3D) optical lattices have been reported PHI ED])- After the submission of this article, 
the observation of the matter-wave DS was reported in ID pT] . 

In contrast to free space, stable localized modes are possible in periodic potentials in 
any dimension both for attractive and repulsive interaction. In the case of a self-repulsive 
BEC generation of multidimensional matter-wave DS due to a modulational instability 
has been predicted theoretically |221; and the existence and stability of 2D DS have 
been studied |231- Using a variational approximation and direct numerical simulation 
the authors of Ref. [21] demonstrated that in the case of attractive interaction above 
the threshold number of atoms, the initial BEC wave packet placed in optical lattice 
collapses into multidimensional DS. 

The effect of the lattice for quantum wave packets much larger than the unit cell 
of the periodic potential, may be replaced by the effective mass. In which case, even 
for repulsive interaction the wave function envelope dynamics, for points in momentum 
space where effective mass is negative, are governed by the NLSE with negative (self- 
focusing) nonlinearity (231201 • The center of the quantum wave packet in momentum 
space can be easily shifted in a controlled manner by accelerating the lattice (for instance 
by chirping the relative detuning of the beams creating the lattice) as was demonstrated 
in early experiments with cold atoms [2Z1 and later with BEC |2H1- Recently, effects 
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of negative effective mass have been studied with ®^Rb condensates in one dimensional 
optical lattices [22] • 

In 2D above critical value of the self-focusing nonlinearity wave packet collapses. 
The nonlinearity in case of BEC determined by number of atoms, scattering length, 
effective mass in the lattice and frequency of transverse confinement. When wave packet 
size becomes comparable to the size of a unit cell effective mass approximation breaks 
down. If the nonlinearity is close to that one for which DS is supported by the band 
gap, part of the wave function is transferred to DS and part decays into linear waves. 
This is a general phenomenon for NLSE when a state is prepared sufficiently close to 
the localized state [HH] . 

In this article, we first obtain a criteria for critical value of interaction using 
variational approximation for the wave packet envelope dynamics |3T| l32] and effective 
mass approximation |2S1 12E] • We also show that this critical nonlinearity is associated 
with the smallest nonlinearity for which the lattice supports DS [331 1211 123, and its 
value is equal to the only nonlinearity with which a stationary solution is possible for the 
free space 2D NLSE with the corresponding effective mass. In the following section, we 
illustrate DS generation with a self-repulsive BEC using two numerical examples: square 
optical lattice for the parameters considered in [221 asymmetric honeycomb lattice. 
The latter was originally considered as a system to study effects of Berry curvature in 
periodic potentials [201 • -^^^ both systems, we simulate all the stages of the possible 
experiments: adiabatic introduction of the lattice, half of the Bloch oscillation and wait 
period for the wave function to collapse to DS. Finally, we summarize the results of this 
study. 

2. Variational Approximation 

The variational approximation for the NLSE was originally developed in nonlinear optics 
(for the review see [221 )• It was successfully applied to describe BEC dynamics [231211 
I2ni mH m] , including its evolution in the optical lattices [HI 1211 112] • In this section we 
first apply variational approximation ideas to the evolution of the BEC Gaussian wave 
packet in free space in D dimensions. 

The Gross-Pitaevskii (CP) equation is a NLSE describing the dynamics of BEC: 



with a nonlinear term being due to the mean-field treatment of the interaction between 
the atoms. It can be both positive and negative depending on the scattering length of 
the atomic collisions. Most of the current experiments deal with self-repulsive BECs 
(positive scattering length). The wave function in the equation is normalized to unity. 

To apply variational approximation, we restrict the dynamics of the quantum wave 
packet in free space (^(r) = 0) in D dimensions to the form 



in—iij{v,t) = Hip{v,t) 



( 



2m 



V^ + V{T) + NgDmT,t)f)^{r,t), (1) 




(2) 



Dynamical generation of two-dimensional matter-wave discrete solitons 



4 



where a and (3 are variational parameters, a being inversely proportional to the width of 
the wave packet squared. The semi-classical Lagrangian corresponding to GP equation 
may be written as 

L(a, d, /3, {3) = {^!d\ th^^ -H\'^d). (3) 
The terms necessary for its calculation are 

(*«|.Sajl'p.> = --, (4) 

The first term describes temporal evolution, the second is associated with kinetic energy, 
and last term describes mean field interaction between particles. After the substitution 
7 = 1/a the Lagrangian equation becomes 

77 - o7 = — T + ——7 Urr • (7) 



In the free space case, without mean field interaction (gf^ = 0), this equation gives an 
exact result for dependence of the wave packet dispersion on time 

f = 2(ao^ + At2), (8) 

where cTq = 7^(t = 0)/2 is initial spatial second moment of the wave packet, and 

~ 2mV(^ = 0) ~ 4m Vg' ^ ^ 

Equation ((Tj) is particularly simple in 2D. In this case, the right hand side becomes 
independent of parameters of the wave packet, and the dynamics is effectively described 
in the same way as for non-interacting atoms. The wave packet becomes dispersionless 
when the right hand side of ((Tj) vanishes. This gives a criteria for critical interaction 
strength 

Ng, = . 10 

m 

For negative nonlinearity, when \g\ > gc, the wave packet collapses. 

Quantum motion of a wave packet in periodic potential can be effectively described 
as motion in free space using the concept of effective mass, which in principle may be 
negative. We comment on this in the next section. 
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3. Effective Mass 

When the external potential V{y) is periodic, and the size of the quantum wave packet 
is much larger than its unit cell, the effect of the potential for different wave vectors ko 
may be described in terms of the effective mass. This mass may be inferred from the 
band dispersion [23121! 

.2/ d^E V' 

When the wave packet is delocalized in real space over many lattice sites, in momentum 
space it is localized around a given wave vector ko. In this case it may be represented 
as a sum of envelopes /n(i") over different Bloch bands 

^(r,t) = 5:/„(r,t)0„ko(r)e-^^-o*/^ (12) 

n 

where the envelope fn{T,t) is a slowly varying function within the unit cell, and each 
Bloch function 0„ko is the solution of the linear eigenproblem 

("I^V^ + V^(r)) 0„k(r) = i?nk0(r), (13) 

normalized to the of the unit cell 

rfr|0„k„(r)|' = a (14) 



cell 



In the experiments it is possible to prepare wave packets that populate only the lowest 
band j2Zl. The NLSE for the envelope incorporates the effects of the external potential 
into effective mass |2S1 IIEI 



« (f + V, . V/„) ^ {-^^^^ + m I//) U, (15) 
where is drift velocity of the wave packet center 

= (</'nko|p|</'nko) , (16) 

and effective interaction strength g2 is given by 

92 = ^^1 rfr|0„ko(r)r- (17) 

cell 

In the particular examples we will discuss below, the wave packet will be driven to the 
point in the Brillouin zone where = 0, and effective mass tensor is negative in all 
directions, hence the envelope dynamics will be governed by NLSE with negative mass. 
This can be viewed as a NLSE with positive mass with inverted sign of nonlinearity, 
which is clearly seen from the equation for complex conjugate of the envelope function 

dt ( 2\meS,fMu\dXf,dx^ ^92\fn\^ fn- (18) 

This equation can be obtained by taking the complex conjugate of (fT3|) and using 
absolute value of the mass. 
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Hence, in this situation, as long as the condition for the effective mass 
approximation holds, m in the Eq. ((7j) should be replaced by nic^. In the case, when 
an effective mass and nonlinearity have opposite signs and \g2\ > g2c-> where critical 
interaction strength is defined in ()10|). with mass being replaced by effective mass, the 
wave packet collapses. Notice that in general from (jZj) it follows that dispersion of the 
wave packets for any D when the effective mass is negative is described by the equation 
with positive mass and inverted sign of nonlinearity. When the size of the wave packet 
becomes comparable to the lattice spacing the effective mass approximation no longer 
holds. As we will discuss in the following section, in 2D there is nonlinearity below 
which the DS are not supported by the bandgap. If the nonlinearity is sufficiently larger 
than this delocalizing nonlinearity, part of the wave packet decays into DS and part 
decays into linear waves. 

4. Delocalization 

In contrast to ID, where DS may correspond to arbitrary nonlinearity, in 2D and higher 
dimensions DS are possible only for nonlinearity above a critical value |HH] • The authors 
of [S3 considered the possibility of observing the delocalization transition with matter- 
wave DS in optical lattices, when an irreversible change from DS to delocalized states 
is produced for a slow change in the lattice parameters. This delocalizing nonlinearity 
may be associated with the critical nonlinearity for a Gaussian wave packet to collapse 
as discussed in Section El 

The concept of the effective mass is not generally applied to the DS since in 
the middle of the gap the DS is localized within one lattice site. As the chemical 
potential of the DS comes close to a band of linear states, its space extension increases, 
hence one may expect that the effective mass approximation becomes applicable. The 
results of Section |21 imply that there is only one value of nonlinearity for which the 
localized modes are supported in 2D free space when the envelopes of the localized 
modes are approximated by Gaussians. This also can be shown in general from the 
scaling arguments for 2D NLSE. Indeed, if the normalized wave function ipi{x,y) is the 
solution of 

- + = /iV^i, (19) 

then the normalized wave function ip2 = BtpilAx, Ay) is the solution of 

-V^tlj2 + 1^'2 = A^jl^Pi. (20) 

This means that in 2D free space there is only one possible value of 7 for which localized 
modes can be found for any value of fl. 

Localized wave packets with very large extension correspond to the critical 
nonlinearity. As their size is reduced, other corrections due to the lattice also start 
to play a role. The variational approximation gives a clear picture of what happens 
to initially localized states of NLSE in 2D. It predicts the critical value of nonlinearity 
above which evolution of the wave packet width changes character. One may expect that 
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the value of Nqc given in (jiup is close to the exact one. We confirmed this expectation by 
performing direct self-consistent numerical simulations based on the effective potential 
approach suggested in jHS]- Similar to the self-consistent Hartree-Fock approximation, 
one may consider a DS to be a localized state in the effective potential created by itself 

V,^{v) = -Ng'o\i^{v)\\ (21) 

In the numerical simulation, we started with an arbitrary nodeless initial wave function, 
and with the imaginary time evolution, found the ground state of the potential (PT|) for 
the Hamiltonian given by 

Hes = + Kff(r), (22) 

2 |meflf| 

and used it in the next step of iteration. We found that independent of the initial guess 
state above the critical value 

5.850 

Ng'num = n T' 23 

\mes\ 

the self-consistent procedure resulted in collapsing states with infinite negative energy, 
while for smaller nonlinearities the states expanded, with energy going to zero. This 
value differs from the one for an extended Gaussian wave packet to collapse (fTUj) by 
~ 10 percent. 

As an alternative method we also reduced the 2D equation to a ID ordinary 
differential equation and solved two point boundary value problem {ip'{0) = 0, tploo) = 
0) with the shooting method We found that for arbitrary energy, nodeless solitons 
in free space with arbitrary size are supported only for one value of nonlinearity given 
by the same value as in ()23|1 . 

The direct dynamical simulations discussed in the next section confirm the existance 
of the critical nonlinearity above which a wave packet with a finite size collapses. Also, 
using the imaginary time evolution for different values of chemical potential in the 
gap, we find corresponding nonlinearity for DS in the lattice to exist. The minimum 
nonlinearity found as a result of this calculation is also in a good agreement with ()10|] 
and 



5. Numerical Simulations 



5.1. General Remarks 

By choosing appropriate units of length, L„, and mass, M„ = Matom, the GP equation (H)) 
is reduced to a dimensionless form 



VL{r) + Ng2\ij\^ + Y-v 



(24) 



We also add an external force F that in the case of optical lattices may be created by 
accelerating the lattice, for example, by sweeping the relative frequency of the beams 
creating the lattice. In this case, unit of energy of the problem is Eu = /MuL"^ and 
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the unit of time is given by t^ = fi/E^. When the unit of length is chosen to be equal 
to the inverse wave vector of the light creating the optical lattice = typical 
maximum depths of the optical potentials achievable in the dissipationless regime are 
~ 20. The forces created by accelerating the lattice are limited due to the finite lifetime 
of the excited states, in case of alkali atoms to about 1000, in the experiments [2Zj 
forces on the order of ~ 1 — 10 were used. The dynamics of the BEC is described 
by a 2D equation in the case when the strong confinement in the transverse direction 
"freezes" the wave function in that direction to the harmonic oscillator ground state. 



This happens when oscillator length in that direction 1^ = y^/Matomt^^ becomes smaller 

than the condensate healing length = (Anna) , where n is atom density and a is 
scattering length. When, at the same time, 4 is still larger than 3D scattering length. 



the collisions between the atoms preserve their 3D character, yet the dynamics of the 
BEC in the two other directions are effectively described by 2D GP equation Such 
a regime was recently demonstrated experimentally [IH]. The nonlinear coefficient in 
this case is given by 



The experimental system, therefore, provides great flexibility with which the nonlinear 
coefficient in corresponding NLSE can be controlled. Parameters variable in experiments 
are number of atoms, A^, the frequency of transverse confinement oOz, and the 
scattering length, a, that in principle may be tuned by a magnetic field with Feshbach 
resonances jlHI. For the experimental parameters of jl^I we estimate the nonlinear 
coefficient in present units to be Ng2 ~ 6000. For the ansatz chosen above to be a good 
choice, kinetic energy (jSJ should be larger than interaction energy It means that 
the nonlinearity should be smaller than ~ 27r. Below we consider only the cases when 
this holds. For BEC the effective nonlinearity can be always reduced either by changing 
number of atoms or trapping frequency in transverse direction. 

The stationary states of Eq.()24j) are described by solutions in the form: ip{r,t) = 
0(r) exp(— i/xt), where /i is the chemical potential. Localized states can be found for /i 
in the gaps of linear problem. 

Below, we consider two examples. The first one is based on parameters for which 
existence and stability of the solitons were studied in |2S1- Our consideration extends the 
treatment to suggest a specific approach for DS generation based on the ideas outlined 
in previous sections. The second example deals with the asymmetric honeycomb lattice 
that was studied by us [3^] in the context of observing self-rotation and Berry curvature 
effects for quantum wave packets in asymmetric periodic potentials. There, robust 
spontaneous generation of the DS above a critical interaction strength was observed for 
wave packets left at the corner of the Brillouin zone. Here, we discuss how this effect 
could be naturally explained in terms of the effective mass concept. 




a|. 



a\ <lz < ^ 



(25) 




(26) 
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In both cases we have performed simulations in a form that mimics possible 
experiments. We start with a Gaussian wave packet in free space with a size that 
is much larger than the unit cell of the potential. The process may be divided into three 
stages: (1) adiabatic introduction of the lattice potential, (2) acceleration of the lattice 
for half of the Bloch oscillation, and (3) a wait period for the wave packet to collapse. In 
the first two stages, adiabaticity is crucial. An adiabatic criteria can be estimated based 
on Landau-Zener formalism for tunneling between two levels with an avoided crossing. 
In this case the probability to tunnel is given by [13 EHj 

P = e.p(-g, (27) 

where 6 is the gap between the states and a is the rate of change of the gap. This 
criteria, when a periodic potential with amplitude Vq is introduced in time ty, becomes 

tv 



6H, 



> 1, (2^ 



Vo 

where 5 is the gap between the first and the second bands at A; = 0, since originally 
the width of the wave packet in the momentum space is much smaller than the size of 
the Brillouin zone (note that units with ^ = 1, are used). The maximum "force" of the 
drive may be similarly estimated as 

>1. (29) 



ok 



Here eg is the smallest gap between the first and the second bands, and || is the largest 
slope of the dispersion. An alternative expression may be derived based on the WKB 
approximation [IH] 

'^-^ » 1 (30) 
F 

where the effective mass, meg, can be estimated from the dispersion as well. We made 
sure that the adiabaticity criteria are fulfilled in the simulations discussed below. We 
also checked it numerically by increasing the driving force and observing the break down 
of the adiabaticity. 



5.2. Square Lattice 

As a first example, we discuss the model potential considered in [221 for existence and 
stability of 2D DS in continuous potentials, namely 

Vz,(x,?/) = Vo(sin^x + sin^y). (31) 

Such a potential may be experimentally obtained by overlapping two pairs of counter- 
propagating beams, far detuned from atomic resonance to reduce effects of spontaneous 
emission. As in [221, we use amplitude Vq = 5. 

In Fig. |l(a)[ we show the dispersion of the linear problem (|T3|) for the potential (j3T|) 
along the high symmetry directions. The first band has a minimum at zero momentum 
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(d) 



Figure 1. (a) - Dispersion for square lattice potential (|31|l along the high symmetry 
path. The green dashed line shows the boundary between strongly bound and quasi- 
unbound states |23| . (b) - xx component of the effective mass tensor along ^^.-axis. (c) 
- Dependence of wave packet dispersion on time after the lattice potential is introduced. 
The two sets of curves are for different points in Brillouin zone: (I) - point F, (II) - 
point M. (d) - Dispersive characteristic x from (|39|l . The blue dots are obtained from 
fitting quadratic dispersion to continuous simulation data from panel (c), red line 
is expected behavior of x- 



(r point), and a maximum at the corner of the Brillouin zone (M point). There exist 
directions on the plane for which the maximum value of the potential is smaller than 
Vq. These are orthogonal {x,y) directions, while the absolute maximum of the potential 
is 2Vo. For the chemical potential of the localized states above Vq (shown as a green 
dashed line in Fig. |l(a)| ), the BEC states are quasi-unbound. For values of /x close to this 
boundary an adequate description of DS is not possible within one-band tight-binding 
model. This is because situations when nodes of the solitonic wave functions are located 
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at the potential minima are possible j2S]- 

The xx-component of the effective mass tensor pip is shown in Fig. |l(b)[ Since 
the potential is separable and symmetric, the ?/?/-component has the same dependence 
on ky, and the tensor is diagonal. At the points of global maximum and minimum of 
the dispersion, the effective mass in both directions is the same. At the V point, the 
effective mass has the smallest positive value (meff^r = 5.53), while at the point M it 
is negative and has the smallest absolute value (meg^M = —5.03). Hence the smallest 
nonlinearity is necessary for the wave packet to self-collapse at the M point. The fact 
that in ID, for a sinusoidal potential for wave vectors, fc, larger than half of the largest 
vector in the Brillouin zone, k > kait = 7r/(2a), where a is lattice spacing, allows to give 
a physical explanation of the origin of the Landau instability studied in [SU] in terms 
of tight-binding approximation. For any interaction strength, wave packets composed 
of Bloch waves for k < /ccrit remain wave packets composed of Bloch waves, while for 
k > kcTit, when interaction is large enough, they partially collapse to localized modes. 

To investigate the validity of the effective mass description, we checked in numerical 
simulations the predictions that may be based on the formulas discussed in previous 
sections. From 0, we introduce the quantity characterizing dispersion of the wave 
packet 

which changes linearly with interaction. As it becomes negative, the wave packet 
collapses. Notice that in this expression the effective interaction strength is related to 
continuous interaction strength with ()17j) . which for the case of separable sinusoidal 
potential at least partially may be computed analytically. For the separable 
potential the solution of the stationary eigenproblem (fT^ is separable: y) = 
4>x{x)4>y{y), where each wave function is given by the solution of the corresponding 
equation 



2(9x2 V 2 y 2 ' ^ ' 

Due to Bloch-Floquet theorem, the solution is a product of a periodic Bloch function 
Uk{x) and a plane wave 

V'x,fc.(a;) = e*'=^^nfc(x). (34) 

After the following substitutions 

b = 2{E-f), 

^x{x) = y{z), 
Eq. becomes the Mathieu equation j^I] 

^ + ib-2qcos2z)y = 0, (36) 
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as a result the solution of (|33p are Mathieu functions with characteristic value h and 
parameter q 

^.,kAx) = Ce(6, g, z) = Ce ^2 - y) , - y, ^) • (37) 
Introducing the numerical factor 

\^ 

(a, g, z)dz 1 

(38) 




7r/2 N 

/ Ce'^.{a,q, z)dz 

-7r/2 y 



we get the following dependence of x oii interaction strength, 

X = 1 + Ng2. (39) 

Therefore the critical interaction strength is given by 

2tt 

Ng2,c = rj- (40) 

For the considered example, area of the lattice unit cell is Q = tt^ and the numerical 
factor for the top of the band is lu = 0.398, which makes the critical interaction strength 
Ng2,c = 0.285. 

To test predicted behaviour of the parameter ^, we perform a numerical simulation 
with a continuous potential. We start with a wave packet of size ax,y = 15/^2, while 
the size of the unit cell is vr in each direction. The lattice potential is ramped from 
Vo = to Vo = 5 in = 450, then the lattice is accelerated with force F = 0.01 
in diagonal [11] direction. The acceleration is halted after the wave packet undergoes 
half of the Bloch oscillation, in other words, when its center is located at point M in 
momentum space. After this, the wave packet expands freely in the presence of the 
lattice potential. Ramp time and the force satisfy adiabaticity conditions (PHj) -()30 |) . so 
that the wave packet during evolutions stays in the first band. In addition to expanding 
the wave packet at the top of the band, we have studied the expansion as soon as the 
lattice was introduced without introducing external field, F. In both cases the agreement 
with analytic prediction (jH^ is very good (see Fig.|T(^and Fig.[T(d)l). 

As the external force is removed the wave packet contracts due to phases 
accumulated during the acceleration even in linear case. But this is only when the 
interaction is above critical that the localized modes are formed. To find a stationary 
localized solution we followed optimization procedure based on a descent technique 
with Sobolev preconditioning used in [22] and described in [22]. In Fig. |2(a)] we display 
dependence of the nonlinearity on chemical potential of the DS. The error bar shows the 
estimated uncertainty for the curve to intersect the band of extended states obtained 
by interpolating to infinite size and infinite relaxation time of the descent procedure. 
The agreement with the argument based on effective mass and free space solitons given 
in Section m is excellent. As one starts with an extended wave packet and nonlinearity 
supported by the gap, it shrinks and may lose some part to radiation (extended states). 
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Figure 2. (a) - Dependence of nonlinearity on corresponding chemical potential. 
The first and second bands are shown with solid rectangles. The inset shows a region 
denoted by a rectangle in the lower left corner. Horizontal solid line shows numerical 
value for delocalizing nonlinearity, gnum, dashed line is the critical nonlinearity for 
extended Gaussian wave packet to collapse 52, c- (b) - Spatial distribution of DS 
corresponding to /i = 0.7 obtained with the descent method, (c) - Probability 
distribution for BEG wave fimction evolved with Ng2 = 0.2 for At ^ 1500 after it 
was driven to M point, (d) - The same for Ng = 0.4 (point A in (a)), approximately 
0.72 of the wave function probability is transferred to the soliton, which corresponds 
to an effective nonlinearity of Ng2 eS ~ 0.288 - at the top of the first band (point B in 
(a)). 



so that the effective interaction experienced by the locahzed mode is then just a fraction 
of the actual interaction. As shown in Fig. |21 the sohtons are formed inside the gap. 
When nonhnearity increases, their chemical potential increases but not significantly, so 
that they stay relatively close to the top of the first band. 
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(c) (d) 

Figure 3. (a) Dispersion for asymmetric lioneycomb potential H41() along high 
symmetry path. The green dashed line shows boundary between strongly bound and 
quasi-unbound states )28| . (b) - xx, and yy components of effective mass tensor along 
fcy-axis. (c) - Dependence of wave packets dispersion on time after the lattice potential 
is introduced. The two sets of curves are for different points in Brillouin zone: (I) 
- point r, (II) - point K. (d) - Dispersive characteristic x from H39|l . The blue dots 
are obtained from fitting quadratic dispersion (jHJ to continuous simulation data from 
panel (c). red line is expected behavior of x- 

5.3. Asymmetric Honeycomb Lattice 

In our previous work jHHj dedicated to study of quantum wave packet propagation 
in spatially asymmetric optical lattices above critical interaction strength, we have 
observed self collapse of the wave packet into localized modes. The effective mass concept 
provides an explanation of the phenomena and relates the critical interaction strength 
to other parameters of the problem. Here, we describe briefly the model potential. The 



Dynamical generation of two-dimensional matter-wave discrete solitons 



15 




(a) 




(c) 



I 
I 



0.018 

0.016 

0.014 

0.012 

0.01 

0.008 

0.006 

0.004 



I 
I 



0.3 



0.25 



0.2 



0.15 



0.1 




(b) 




(d) 



I 



0.6 



0.5 



0.4 



: 0.3 

1^0.2 



I 



I 



0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 

0.04 



Figure 4. Probability density of the wave packet that has been driven to K point 
of square lattice after expansion for Atoxp = 200 with different interactions: (a) - 
Ng2 = 0.5, (b) - Ng2 = 2.5 (c) - Ng2 = 5 (d) - Ng2 = 10. As the interaction increases 
fraction of the wave function transferred to solitons decreases. 



potential used in [36] is 

6 

Mat(r) = ^ A cos(k, ■ r + (41) 

i=l 

with the following parameters 

= ^3 = ^5 = -4, 

A2 = Ai = Ae = -4.664, ^^2) 

01 = 03 = 05 = 0, 
02 = -04 = 06 = l-l- 

It may be obtained either with interference of three pairs of counter propagating beams 
or by holographic techniques. The band structure is shown in the Fig. |3(a)| Asymmetry 
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splits the first band in the middle. In this case, the boundary between strongly bound 
and quasi-unbound states is located in the second gap. The effective mass is negative 
in both directions at K point (see Fig. |3(b)| ). We start with a Gaussian wave packet 
in free space which has spatial dispersion a^^y = 10/a/2 ^ 7.071. Since the size of the 
unit cell in this case is = 3\/3/2, the wave packet occupies several unit cells. The 
lattice potential is introduced in time ty = 120 and for we accelerate it along y-axis with 
|F| = 0.05. In this case conditions in ()28|) - (jnn|) are fulfilled. Analysis of the expansion 
for different values of the interaction at the bottom of the first band (F point) and 
the top of the first band (K point) is shown in the Fig |3(c)| To compare results of 
the continuous simulation with the prediction of the ()32|1 we calculated effective mass, 
Mgff, from the band structure and the numerical factor /, that involves integrals of 
the Bloch wave functions, by integrating over one unit cell wave functions obtained by 
adiabatic evolution. Their values are Mcs^r = 1.7986, Jr = 0.9396 at the point F and 
Mes,K ~ —0.8918, Ik ~ 1.9567 at the point K. Probability density plots after expansion 
for Atexp = 200 at the K point are shown in the Fig. |3] As interaction increased past 
critical interaction the lattice solitons are dynamically formed. The figure illustrates that 
the phenomena is also observed for a square lattice: the fraction of the wave function 
transferred to the localized modes decreases as the interaction strength is increased. 
The effective interaction experienced by each mode is such that corresponding chemical 
potential is close to the first band. 

As we have seen in the numerical examples discussed above, if the nonlinearity is 
small enough, the wave packet collapses into a single DS. Clouds of ultracold atoms can 
be imaged non-destructively [22] • Observation of a persistent atomic cloud localized 
to dimensions comparable to the wave length of the light forming the lattice and will 
be a clear signature of the DS formation. Observations of the predicted delocalizing 
transitions with a single DS prepared in a controlled fashion for a varying lattice depth 
is also an exciting possibility. 

6. Conclusion 

In this article we have shown how DS could be generated with a self- repulsive BEG 
in optical lattices by driving the wave packets to the points where the effective mass 
is negative, that leads to self collapse of the wave packet. As a result, part of the 
wave function is transferred to localized modes with chemical potential in between the 
first and second bands where DS with positive nonlinearity are supported. In 2D, this 
happens only when the critical nonlinearity is achieved. As the strength of interaction 
increases, the wave packet may excite several localized modes. We have observed that 
effective interaction of each mode is such that the corresponding chemical potential is 
located close to the top of the first band. 
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